# -*- coding: utf-8 -*-
"""
Analyse breakjunction data
"""
from mayavi import mlab
from os.path import isfile
import scipy
import scipy.linalg
import pylab
import pickle
from threading import Thread
from time import sleep

from molel.tools.calc import rollavg, reduceData, mycorrcoef
from molel.tools.progressbar import progressbar 
from molel.tools.plotroutines import plotIMs, my_cmap, plotim, plotims, mysavefig
from molel.loaddatfile import loaddatfile
import gzip

class Slopediff(Thread):
    """
    What am I doing
    """
    def __init__ (self, matrix_ai, matrix_a, lnbreake):
        Thread.__init__(self)
        self.matrix_ai = matrix_ai
        self.matrix_a = matrix_a
        self.lnbreake = lnbreake
        self.residuals = None
        self.par = None
    def run(self):
        self.par = self.matrix_ai*scipy.matrix(self.lnbreake)
        self.residuals = (scipy.array(self.matrix_a*self.par-scipy.matrix(self.lnbreake))**2).sum(0)


class FindBreakPoints():
    """
    Finds the rupure of the last breakpoint in a breaking data.
    """
    def __init__(self, g_min_break=0.9, d_g_threashold=.4, g_threashold=1.):
        self.g_min_break = g_min_break
        self.d_g_threashold = d_g_threashold
        self.g_threashold = g_threashold
    def __call__(self, breake, make):
        stepstart_breake = pylab.find(breake[:]>self.g_threashold)
        if len(stepstart_breake)>0:
            stepstart_breake = stepstart_breake[-1]
        else:
            stepstart_breake = 0
        i = None 
        for i in xrange(stepstart_breake, len(breake)):
            try:
                d_g_tmp = breake[i+1]/breake[i]
            except IndexError:
                break
            if (d_g_tmp < self.d_g_threashold)    or (breake[i] < self.g_min_break):
                break
        stepstart_breake = i
        stepstart_make = pylab.find(make>self.g_threashold)
        if len(stepstart_make)>0:
            stepstart_make = stepstart_make[0]
        else:
            stepstart_make = 0
        return stepstart_breake, stepstart_make

        
class Bjdata(object):
    def __init__(self, **kw):
        """
        MCBJ data
        """
        self.breake = scipy.nan
        self.motor = scipy.nan
        self.make = scipy.nan
        self.meta = scipy.nan
        self.delta_opening = scipy.nan
        self.filenames = scipy.nan
    def savez(self, filename):
        """
        Save data, and results to filename
        """
        pickle.dump(self.__dict__, gzip.open(filename, 'wb',5))
        #scipy.savez(filename, breake=self.breake, make=self.make, meta=self.meta,
        #            motor=self.motor,delta_opening=self.delta_opening)
    def loadz(self, filename):
        """
        Load data, and results from filename
        """
        
        self.__dict__ = pickle.load(gzip.open(filename, 'rb'))
        #tmp = scipy.load(filename)
        #map(lambda x: setattr(self, x, tmp[x]), tmp.files)
        #self.meta = self.meta.tolist()
        
    def loadfiles(self, filenames, resizedata=1., autoresize=True, 
                  force=False, **kw):
        """
        Load experimental files (one per breaking/making curve) which matches filenames. 
        If force is true and a "dat.gz" file has not been generated it loads from the 'dat'
        files, otherwise it loads the 'dat.gz' file
        """
        self.filenames = filenames
        if isfile(filenames.replace('*.dat', '.dat.gz')) and (not force):
            self.loadz(filenames.replace('*.dat', '.dat.gz'))
            print ("Import from %s. To import from dat files use force=True"
                   % filenames.replace('*.dat', '.dat.gz'))
            return True
            
        import commands
        filer = commands.getstatusoutput("ls -t -r "+str(filenames))[1].split("\n")
        largest_filer = commands.getstatusoutput("ls -S "+str(filenames))[1].split("\n")[1]
        self.motor = []
        self.meta, rdat = loaddatfile(largest_filer, **kw)
        kw = self.beforeimport(kw)
        self.breake = []
        self.make = []
        if(autoresize):
            resizedata = round(len(rdat)/5000.)  
        
        pbar = progressbar(len(filer))
        for fil in filer:     #Load s All files containing breaking / make traces into list 'data' 
            pbar()
            meta, rdat = loaddatfile(fil, **kw)
            if rdat.ndim == 2 and rdat.shape[0]>50:
                if resizedata > 1:
                    rdat = reduceData(rdat, resizedata)
                self.make.append(scipy.array(filter(lambda a: a != 0, rdat[:, 2])))
                self.breake.append(scipy.array(filter(lambda a: a != 0, rdat[:, 1])))
                self.motor.append(meta['Motor_pos[nm]'])
                self.delta_opening = (rdat[1, 0]-rdat[0, 0])   # in opening
        self.motor = scipy.array(self.motor)
        self.correctdata()
        self.savez(self.filenames.replace('*.dat', '.dat.gz')) 
    def beforeimport(self,kw):
        """
        
        """
        return kw        
    def correctdata(self):
        """
        Inplement correction of data (e.g. correction for extra resistance)
        """
        return None    
       
class Bjanalyse(object):
    def __init__(self, bjdata, **kw):
        """
        MCBJ analysis
        """
        self.bjdata = bjdata
        self.g_bins = []
        self.d_length = None
        self.his_breake = scipy.nan
        self.his_make = scipy.nan
        scipy.seterr(divide='ignore')
        self.sub = None
        
    def savez(self, filename):
        """
        Save data, and results to filename
        """
        pickle.dump(self.__dict__, gzip.open(filename, 'wb',5))
        
    def loadz(self, filename):
        """
        Load data, and results from filename
        """
        
        self.__dict__ = pickle.load(gzip.open(filename, 'rb'))
   
    def findbreakpoints(self, findbreakpoint_func=None, **kw):
        self.breakpoints, self.makepoints = [], []
        findbreakpoint_func = findbreakpoint_func if (findbreakpoint_func is not None) else FindBreakPoints()
        for i in range(len(self.bjdata.breake)):
            (breakpoints_t, makepoints_t) = findbreakpoint_func(self.bjdata.breake[i], self.bjdata.make[i], **kw)
            self.breakpoints.append(int(breakpoints_t))
            self.makepoints.append(int(makepoints_t))
        self.breakpoints = scipy.array(self.breakpoints)
        self.makepoints = scipy.array(self.makepoints)
    def twodimdata(self):
        max_breakepoint = self.breakpoints.max()
        max_makepoint = self.makepoints.max()
        datalength_breake = (scipy.array(map(len, self.bjdata.breake))-scipy.array(self.bjdata.breakpoints)).max()
        datalength_make = (scipy.array(map(len, self.bjdata.make))-scipy.array(self.bjdata.makepoints)).max()
        twodim_breake = scipy.zeros((datalength_breake+max_breakepoint, len(self.bjdata.breake)))
        twodim_make = scipy.zeros((datalength_make+max_makepoint, len(self.bjdata.make)))
        for i, p in enumerate(self.breakpoints):
            twodim_breake[(max_breakepoint-p):(max_breakepoint-p+len(self.bjdata.breake[i])), i] = self.bjdata.breake[i]
        for i, p in enumerate(self.makepoints):
            twodim_make[(max_makepoint-p):(max_makepoint-p+len(self.bjdata.make[i])), i] = self.bjdata.make[i]
        return twodim_breake, max_breakepoint, twodim_make, max_makepoint
        
    def cycledata(self):
        self.cycle_breakedata = [[], []]
        self.cycle_makedata = [[], []]
        for i in range(len(self.bjdata.make)):
            self.cycle_breakedata[0] += list(self.bjdata.motor[i,1] - scipy.arange(len(self.bjdata.breake[i]))*self.bjdata.delta_opening)
            self.cycle_makedata[0] += list(self.bjdata.motor[i,0] + scipy.arange(len(self.bjdata.make[i]))*self.bjdata.delta_opening)
            self.cycle_breakedata[1] += list(self.bjdata.breake[i])
            self.cycle_makedata[1] += list(self.bjdata.make[i])
                
    def plot_cycle(self):
        if not self.__dict__.has_key('cycle_breakedata'):
           self.cycledata()
        bins = (scipy.linspace(min(self.cycle_breakedata[0]), max(self.cycle_breakedata[0]), 200), 
               10**scipy.linspace(scipy.real(scipy.log10(min(self.cycle_breakedata[1]))), 
                                  scipy.real(scipy.log10(max(self.cycle_breakedata[1]))), 200))
        hist = scipy.histogram2d(self.cycle_breakedata[0], self.cycle_breakedata[1], bins=bins)
        pylab.subplot(121)
        pylab.imshow(scipy.log10(hist[0].T), origin='lower', 
                     extent = [min(self.cycle_breakedata[0]), max(self.cycle_breakedata[0]), 
                               scipy.real(scipy.log10(min(self.cycle_breakedata[1]))),scipy.real(scipy.log10(max(self.cycle_breakedata[1])))], 
                     aspect='auto')
        pylab.xlabel('um')
        pylab.title('Breaking')
        bins = (scipy.linspace(min(self.cycle_makedata[0]), max(self.cycle_makedata[0]), 200),
               10**scipy.linspace(scipy.real(scipy.log10(min(self.cycle_makedata[1]))), 
                                  scipy.real(scipy.log10(max(self.cycle_makedata[1]))), 200))
        hist = scipy.histogram2d(self.cycle_makedata[0], self.cycle_makedata[1], bins=bins)
        pylab.subplot(122);pylab.imshow(scipy.log10(hist[0].T), origin='lower', 
                                        extent = [min(self.cycle_makedata[0]), max(self.cycle_makedata[0]), 
                                                 scipy.real(scipy.log10(min(self.cycle_makedata[1]))), scipy.real(scipy.log10(max(self.cycle_makedata[1])))], 
                                        aspect='auto');pylab.xlabel('um');pylab.title('fusing')

    def imshow(self, breake=True, slog=False, trange = [0, -1], from_g0=False, **kw):
        if not kw.has_key('vmin'):
            kw['vmin'] = -6 if slog else 1e-6         
        if not kw.has_key('vmax'):
            kw['vmax'] = 1 if slog else 1e1        
        if not kw.has_key('interpolation'):
            kw['interpolation'] = 'nearest'
        twodim_breake, max_breakepoint, twodim_make, max_makepoint = self.twodimdata()
        if breake:
            if from_g0:
                trange[0] = max_breakepoint
            y = twodim_breake[trange[0]:trange[1], :]
        else:
            if from_g0:
                trange[0] = max_makepoint
            y = twodim_make[trange[0]:trange[1], :]
        y = y[::(scipy.shape(y)[0]/1000+1), ::(scipy.shape(y)[1]/1000+1)]
        if slog:
            y = scipy.log10(y)
            y[ y == -scipy.Inf ] = -10
        pylab.imshow(scipy.real(y), **kw)
        
    def correlate_metallicwire(self, **kw):
        makbreak= max(map(scipy.nanmax,self.bjdata.breake))
        self.correlate(scipy.linspace(0., 20. if makbreak>20.  else makbreak, 150),
                     norm_below=10, scale_clim=0.5, **kw)
        
    def select_correlate(self, I, J, threshold=3., nmean=5):
        gtilde = (self.hist_alltraces.T-self.hist_alltraces.mean(1)).T
        tmp = scipy.array([((gtilde*gtilde[i]).T/(abs(gtilde).mean(1)*abs(gtilde[i]).mean())).T for i in range(len(gtilde))]).transpose((2, 0, 1))
        datah = tmp[:, I-nmean//2:I-nmean//2+nmean, J-nmean//2:J-nmean//2+nmean].mean(2).mean(1)
        pylab.clf();pylab.hist(datah, bins=50)
        pylab.figure();pylab.clf()
        return pylab.find(threshold<datah)

    def correlate(self, bins, data = None, rollavg_n=1, normtointerval=.95, norm_below=.1, scale_clim=1.):
        data = self.bjdata.breake if (data is None) else data
        self.g_bins = bins
        codict = {'red': ((0.0, 0.0, 0.0), 
                                    (0.2, 0.2, 0.2), 
                                 (0.32, 0.0, 0.0), 
                                 (.4, 0.8, 0.8), 
                                 (.6, .7, .7), 
                                 (1, 1.0, 1.0)), 
                 'green': ((0.0, 0.0, 0.0), 
                                     (.32, 1.0, 1.0), 
                                     (.4, 0.8, 0.8), 
                                     (.6, 0.5, 0.5), 
                                     (.8, 0.0, 0.0), 
                                     (1, 0.0, 0.0)), 
                 'blue': ((0.0, .4, .4), 
                                    (0.25, 1.0, 1.0), 
                                 (0.32, 0.0, 0.0), 
                                 (0.6, 0.0, 0.0), 
                                 (.8, 0.0, 0.0), 
                                 (1, 0.0, 0.0))}
        islog = scipy.any(scipy.diff(scipy.diff(bins))>1e-14)
        logfunc = scipy.log10 if islog else lambda x:x
        coor_cmap = pylab.matplotlib.colors.LinearSegmentedColormap('my_colormap', codict, 256)
        self.hist_alltraces = scipy.array(map(lambda x: scipy.histogram(x[scipy.isnan(x)==False], bins=bins)[0], data)).T
        g_data = self.hist_alltraces if rollavg_n == 1 else rollavg(self.hist_alltraces, rollavg_n, axis=0)
        corr_data = scipy.corrcoef(g_data)
        corr_data2 = scipy.cov(g_data)
        filt = scipy.tri(len(corr_data), len(corr_data), 0, dtype=int)==1
        test = corr_data2[bins[:-1] < norm_below, bins[:-1] < norm_below].ravel()
        normalize = scipy.sort(test)[int(len(test)*normtointerval)]
        corr_data2 *= 1./normalize
        corr_data[filt] = corr_data2[filt]
        self.corr_data = corr_data
        extent = (lambda x: (max(x), min(x), min(x), max(x)))(logfunc(bins) if max(scipy.diff(bins, 2))>1e-8 else bins)
        pylab.imshow(corr_data[::-1, ::-1].T, origin='upper', extent=extent, vmin=-.5*scale_clim, vmax=1.*scale_clim, cmap=coor_cmap, interpolation='nearest')
        cbar = pylab.colorbar(shrink=.9)
        ax = pylab.gca()
        if islog:
            ax.yaxis.set_ticklabels (map(lambda x: r'$\mathdefault{10^{%i}}$' % x, ax.yaxis.get_ticklocs()))
            ax.xaxis.set_ticklabels (map(lambda x: r'$\mathdefault{10^{%i}}$' % x, ax.xaxis.get_ticklocs()))
        pylab.xlabel("Conductance $G/G_0$")
        pylab.ylabel("Conductance $G/G_0$")
        ax.xaxis.set_label_position('top')
        ax.xaxis.set_ticks_position('top')
        cbar.ax.text(0, 1.02, r'($\times%.1f$)'%normalize, transform = cbar.ax.transAxes, va='bottom', ha='left');pylab.draw()
        
    def platolength_im(self, plotodef_conddecade=.5, resolution=100., **kw):
        """
        Plots plato length vs. conductance.
        plotodef_conddecade: Definition of plato is a datapoints which fall into a conduction window of [10**C; 10**(C-plotodef_conddecade)].
        resolution: Appr. 'x' resolution of picture. Maximum is automatic numbers of datapoints in the longest plato.
        """
        j = int(round(plotodef_conddecade/(scipy.log10(self.g_bins[1])-scipy.log10(self.g_bins[0]))))
        tmp = self.hist_alltraces.cumsum(0)
        new_his = tmp[j:]-tmp[0:-j]
        bins_l = scipy.arange(new_his.min(), new_his[self.g_bins[:-j-1] < 1].max(), 1)
        bins_l = bins_l[:: (round(len(bins_l)/resolution) if round(len(bins_l)/resolution)>0 else 1)]
        im = scipy.array(map(lambda x: scipy.histogram(x, bins=bins_l)[0], new_his))
        bins_l = self.bjdata.delta_opening*bins_l
        pylab.imshow(scipy.log10(im), origin='lower', extent= (bins_l[0], bins_l[-1], scipy.log10(self.g_bins[0]), scipy.log10(self.g_bins[-1])), aspect='auto', **kw)
        pylab.xlabel("Plato length [nm]");pylab.ylabel("Conductance $\log(G/G_0)$")
    def findcut(self, b, dcut=5e-5, dirr=-1):
        """
        return cut in dataseries (to filter length traces on the length 'b' at a C dcut)
        """
        try:
            i = (b>dcut).tolist()[::dirr].index(True)
        except ValueError:
            i = 0
        return i 
    def make_subs(self,  breake=True, pars=(2, 1e-3), force=False):
        """pars = (platoplength, conductance)"""
        if (not force) and ((breake and (vars(self).has_key('his_breakesub'))) or ((not breake) and (vars(self).has_key('his_makesub')))):
            return 0
        if breake:
            n = scipy.array(map(lambda x: self.findcut(x, dcut=pars[1]) , self.bjdata.breake ))-self.breakpoints
            sub = [pylab.find(n <= (pars[0]/(self.bjdata.delta_opening)))]
            sub.append(pylab.find(n > (pars[0]/(self.bjdata.delta_opening))))
        else:
            n = scipy.array(map(lambda x: self.findcut(x, dcut=pars[1], dirr=1) , self.bjdata.make))-self.makepoints
            sub = [pylab.find(n <= -(pars[0]/(self.bjdata.delta_opening)))]
            sub.append(pylab.find(n > -(pars[0]/(self.bjdata.delta_opening))))
        for i in range(len(sub), 0, -1):
            if len(sub[i-1]) == 0: 
                del(sub[i-1])
        self.sub = sub   
        
    def _make_histograms(self, sub=None, breake=True):
        data = self.bjdata.breake if breake else self.bjdata.make
        sub = range(len(data)) if sub is None else sub
        zerolengthpoint = self.breakpoints if breake else self.makepoints
        conducdata = scipy.concatenate([data[i] for i in sub])
        lengthdata = scipy.concatenate([scipy.arange(len(data[i])) - zerolengthpoint[i] for i in sub])
        notnans = scipy.isnan(conducdata*lengthdata)==False
        conducdata = conducdata[notnans]
        lengthdata = lengthdata[notnans]
        self.length_bins = scipy.arange(lengthdata.min(),lengthdata.max(),self.d_length)
        return scipy.histogram2d(conducdata,lengthdata, bins=[self.g_bins,self.length_bins])[0].T*self.d_length*self.bjdata.delta_opening/float(len(sub))
        # normalized to number traces and to the inverse size of lengthbins

    def make_histograms(self,g_bins, d_length=0.01, force=False,  breake=True):
        self.g_bins = g_bins
        d_length = d_length//self.bjdata.delta_opening+1
        if ( force or (self.d_length != d_length) or
             (breake and (scipy.shape(self.his_breakesub[0])[1]!=(len(self.g_bins)-1) if vars(self).has_key('his_breakesub') else True)) or
             ((not breake) and (scipy.shape(self.his_makesub[0])[1]!=(len(self.g_bins)-1) if vars(self).has_key('his_makesub') else True))
           ):
            self.d_length = d_length
            his_sub = []
            his_sublen = []
            if self.sub is not None:
                for s in self.sub:
                    his_sub.append(self._make_histograms(sub=s, breake=breake))
                    his_sublen.append(len(s))
            else:
                his_sub.append(self._make_histograms(sub=None, breake=breake))
                his_sublen.append(len(self.bjdata.breake))
            his = scipy.array(map(lambda x, y:x*y, his_sub, his_sublen)).sum(0)/scipy.sum(his_sublen)
            if breake:
                self.his_breakesub = his_sub
                self.his_breakesublen = his_sublen
                self.his_breake = his
            else:
                self.his_makesub = his_sub
                self.his_makesublen = his_sublen
                self.his_make = his
    def CalcLengthPlateaus(self, Glim = [.5, 1.2], bins=100):
        self.plateaus_lengths = scipy.array(map(lambda data: scipy.sum((Glim[0]<data)*(Glim[1]>data), axis=0),self.bjdata.breake))*self.bjdata.delta_opening*10
        pylab.hist(self.plateaus_lengths, bins=bins, color='k')
        #pylab.text(0.6, .8, (r"$%.2e}>G>%.2e}$" % tuple(Glim)).replace('e+00', r'{').replace('e-00', r'{').replace('e-0', r'\cdot 10^{-').replace('e+0', r'\cdot 10^{+').replace('e+', r'\cdot 10^{+').replace('e-', r'\cdot 10^{-'), transform=pylab.gca().transAxes)
        pylab.title((r"$%.2e}>G>%.2e}$" % tuple(Glim)).replace('e+00', r'{').replace('e-00', r'{').replace('e-0', r'\cdot 10^{-').replace('e+0', r'\cdot 10^{+').replace('e+', r'\cdot 10^{+').replace('e-', r'\cdot 10^{-'))
        pylab.xlabel('Aprox. plateau length [$\AA$]');pylab.ylabel('Trace counts')
    def SnapVsPlateaulength(self, c_lim=(.5, 1.2), snapthreashold=.2, n=20):
        plateaus_lengths = scipy.array(map(lambda data: scipy.sum((c_lim[0]<data)*(c_lim[1]>data))*self.bjdata.delta_opening*10, self.bjdata.breake))
        snap = scipy.array(map(lambda x,x0: x[x0:][pylab.find(x[x0:]<snapthreashold)[0]], self.bjdata.breake,self.breakpoints))
        hist_snap = scipy.histogram(snap, bins=10**scipy.linspace(-6, .1, n))
        list_of_plateaus_lengths = [plateaus_lengths[(snap<(hist_snap[1][i+1]))*(snap>s)] for i, s in enumerate(hist_snap[1][:-1])]
        pylab.plot(hist_snap[1][:-1], map(scipy.median, list_of_plateaus_lengths), 'b', lw=2, drawstyle='steps-post')
        pylab.plot(hist_snap[1][:-1], map(scipy.mean, list_of_plateaus_lengths), 'r', lw=2, drawstyle='steps-post')
        pylab.plot(snap, plateaus_lengths, '.k');pylab.axis(xmin=(10**-6 if pylab.xlim()[0]<=0 else pylab.xlim()[0]));pylab.xscale('log')
        pylab.ylabel('Plateau Length [$\\AA$]');pylab.xlabel('Snap to G');pylab.xlim(None, .1)    

    
    def slope_diff(self, lnm=0.1): 
        Dt = self.bjdata.delta_opening
        data = self.bjdata.breake
        nfitlength = int(lnm/Dt)
        A = scipy.ones((nfitlength, 2))
        A[:, 0] = scipy.arange(nfitlength)*Dt
        A = scipy.matrix(A)
        try:
            Ai = scipy.linalg.inv(A.T*A)*A.T
        except scipy.linalg.LinAlgError:
            print "LinAlgError: singular matrix"
            return False
        maxlendata = max(map(len,data))
        diffmatrix = scipy.zeros((maxlendata,maxlendata))
        Avector = scipy.concatenate((Ai[0,:],scipy.zeros((1,maxlendata-nfitlength))),axis=1)
        for i in xrange(maxlendata):
            diffmatrix[i,:] = scipy.roll(Avector,i-nfitlength/2)
        dI=[]
        for dat in data:
            lnbreake = scipy.log10(dat).astype('float')
            lnbreake[ lnbreake == -scipy.Inf ] = scipy.nan
            lnbreake = scipy.matrix(lnbreake)
            dI.append(scipy.ravel(diffmatrix[:lnbreake.shape[1],:lnbreake.shape[1]]*lnbreake.T))
            dI[-1][range(nfitlength)+range(-nfitlength,0)] = scipy.nan
        return dI 
    def plot_slopes(self, g_bins_lim = (-6, 1), binslen=500, lnm=0.1, d_c_lim=(-15, 5),constrict_d_c=(0, 5), **kw):
        """
        Test
        """
        dI = self.slope_diff(lnm=lnm)
        I = self.bjdata.breake
        dI = scipy.concatenate(dI)
        I = scipy.real(scipy.log10(scipy.concatenate(I)))
        notnans = scipy.isnan(I*dI)==False
        hist_i_vs_di = scipy.histogram2d(I,dI,bins = [scipy.linspace(g_bins_lim[0], g_bins_lim[1], binslen),
                                                      scipy.linspace(d_c_lim[0], d_c_lim[1], binslen)])
        f0 = pylab.figure( figsize=(9, 6),dpi=150)    
        ax1 = pylab.axes([.1,.1,.55,.8])
        if scipy.nansum(hist_i_vs_di[0])!=0:
            pylab.imshow(scipy.log10(hist_i_vs_di[0][:, :]),origin='lower', extent = [d_c_lim[0], d_c_lim[1], g_bins_lim[0], g_bins_lim[1]], 
                         aspect='auto', **kw)
        [pylab.plot([constrict_d_c[i]]*2,g_bins_lim,'k:',lw=2) for i in range(2)]
        pylab.xlabel("$dC/dl$")
        pylab.ylabel("$\log(G/G_0)$")
        constrict_to_bins = (hist_i_vs_di[2]>constrict_d_c[0])*(hist_i_vs_di[2]<constrict_d_c[1])
        hist_i_vs_di[0][:,constrict_to_bins[1:]].sum(1)
        pylab.axes([.65,.1,.3,.8], sharey=ax1)
        pylab.plot(hist_i_vs_di[0][:,constrict_to_bins[1:]].sum(1),hist_i_vs_di[1][1:] )
        pylab.xlabel(r"Counts")
        pylab.ylabel("$\log(G/G_0)$")
        
    def plot2Dhist(self, force=False,  xval=None, yval=None, breake=True, sub=False, log=False, 
                   plot_colorbar=False, Xlim=None, Plot1D=True, **kw):
        logC = any(scipy.diff(self.g_bins, 2)>1e-14)
        plotD = (self.his_breake if breake else self.his_make).T
        x1, x2 = (0, -1)
        xvalr = [self.length_bins[0]*self.bjdata.delta_opening*10,self.length_bins[-1]*self.bjdata.delta_opening*10]
        if xval is not None:
            dY, dX = scipy.shape(plotD)
            dX = (xvalr[1]-xvalr[0])/float(dX)
            if xval[0] < xvalr[0]:
                xval = (xvalr[0], xval[1])
            if xval[1] > xvalr[1]:
                xval = (xval[0] , xvalr[1]) 
            x1, x2 = (scipy.array(xval)-xvalr[0])/dX
            xvalr = xval
        if yval is None:
            yval = [scipy.log10(self.g_bins.min()), scipy.log10(self.g_bins.max())]    if logC else [self.g_bins.min(), self.g_bins.max()]
        #kw['aspect'] = ((xvalr[1]-xvalr[0])/(yval[1]-yval[0]))/( kw['aspect'] if kw.has_key('aspect') else 1.)
        kw['aspect'] = ((xvalr[1]-xvalr[0])/(yval[1]-yval[0]))/( kw['aspect']) if kw.has_key('aspect') else 'auto'
        axparent = pylab.gca()
        v = axparent.get_position()
        dat = (self.his_breakesub if breake else self.his_makesub) if sub else [plotD.T]
        nn = len(dat)
        dy = (v.ymax-v.ymin)/nn - (nn/2)*0.05
        for i, dd in enumerate(dat):
            im = dd
            his = scipy.nansum(im, 0)
            ax1 = pylab.axes([v.xmin, v.ymin+i*(v.ymax-v.ymin)/nn+(nn/2)*0.05*i, (v.xmax-v.xmin)*(0.7 if Plot1D else 1.), dy])
            im = im.T[:, x1:x2] 
            ax1.imshow(scipy.log10(im) if log else im, origin='lower', extent=(xvalr[0], xvalr[1], yval[0], yval[1]), **kw)
            kw['vmax'] = kw['vmax'] if kw.has_key('vmax') else (scipy.ceil(scipy.log10(scipy.nanmax(im))) if log else scipy.nanmax(im))
            kw['vmin'] = kw['vmin'] if kw.has_key('vmin') else (scipy.floor(scipy.log10(scipy.nanmin(im[im>0]))) if log else 0)
            pylab.xlabel('Approx. gap size / $\AA$')
            pylab.ylabel(r"$\log(G/G_0)$" if logC else r"$G/G_0$")
            if logC:
                ax1.yaxis.set_ticklabels (map(lambda x: r'$\mathdefault{10^{%i}}$' % x, ax1.yaxis.get_ticklocs()))
            ax1.set_visible(True)
            if plot_colorbar:
                if breake:
                    ax2=pylab.axes([v.xmin+(v.xmax-v.xmin)*(0.5 if Plot1D else 0.8), v.ymin+i*(v.ymax-v.ymin)/nn+.1*dy+(nn/2)*0.05*i, (v.xmax-v.xmin)*0.02, dy*.8])
                else:
                    ax2=pylab.axes([v.xmin+(v.xmax-v.xmin)*0.1, v.ymin+i*(v.ymax-v.ymin)/nn+.1*dy+(nn/2)*0.05*i, (v.xmax-v.xmin)*0.02, dy*.8])
                cc=pylab.colorbar(ax=ax1, cax=ax2, format="%.1f" if log else '%.1e' )
                if log:
                    for tick in cc.ax.yaxis.get_ticklabels():
                        tick.set_color('black')
            if Plot1D:
                ax3=pylab.axes([v.xmin+(v.xmax-v.xmin)*0.7, v.ymin+i*(v.ymax-v.ymin)/nn+(nn/2)*0.05*i, (v.xmax-v.xmin)*0.3, dy])
                y=(self.g_bins[:-1]+self.g_bins[1:])/2. 
                if scipy.nansum(his)!=0:
                    ax3.plot(his, y);pylab.yscale('log' if logC else 'linear')
                    pylab.axis(xmin=(.1 if pylab.xlim()[0]<=0 else pylab.xlim()[0]));pylab.xscale('log')
                    if Xlim is None:
                        pylab.xlim(10**scipy.ceil(scipy.log10(scipy.nanmax(his[y<1]))), 10**scipy.nanmin([scipy.floor(scipy.log10(scipy.nanmin(his[y<1]))), scipy.log10(scipy.nanmax(his[y<1]))-1]))
                    else:
                        pylab.xlim(Xlim)
                    ax3.xaxis.set_ticks_position ('top')
                    pylab.ylim((10**yval[0], 10**yval[1]) if logC else    (yval[0], yval[1]))
                ax3.yaxis.set_ticks_position ('right');pylab.grid(True)
            axparent.set_visible(False);ax1.set_visible(True)
    def histogram(self, force=False, bins=10**scipy.arange(-5, 1, 0.02), bars=True, breake=True, scale='log', sub=False, **kw):
        self.make_histograms(g_bins=bins,force=force,  breake=breake)
        width = self.g_bins[1:]-self.g_bins[:-1]
        if sub:
            v = pylab.gca().get_position()
            pylab.gca().set_visible(False)
            dat = self.his_breakesub if breake else self.his_makesub
            nn=len(dat)
            for i, dd in enumerate(dat):
                pylab.axes([v.xmin+i*(v.xmax-v.xmin)/nn, v.ymin, (v.xmax-v.xmin)/nn-0.05, (v.ymax-v.ymin)])
                dat = dd.sum(0)
                pylab.bar(self.g_bins[:-1], dat , width=width, linewidth=0, **kw) if bars else pylab.plot(self.g_bins[:-1]+width/2., dat, **kw)
                pylab.xscale(scale)
            pylab.ylim(0, (pylab.gca().lines[-1].get_ydata ()[self.g_bins[:-1]<1]).max())
        else:
            dat = self.his_breake.sum(0)    if breake else self.his_make.sum(0)
            pylab.bar(self.g_bins[:-1], dat , width=width, linewidth=0, **kw) if bars else pylab.plot(self.g_bins[:-1]+width/2., dat, **kw)
            pylab.xscale(scale)
            pylab.ylim(0, (pylab.gca().lines[-1].get_ydata ()[self.g_bins[:-1]<1]).max())
    def Plot2D_corr(self, binlength=100, normtointerval=.95, rollavg_n=1, scale_clim=1, **kw):
        parentax= pylab.gca()
        v=parentax.get_position()
        pylab.axes([v.xmin+0.12*(v.xmax-v.xmin), v.ymin+0.15*(v.ymax-v.ymin), .25*(v.xmax-v.xmin), .8*(v.ymax-v.ymin)])
        self.plot2Dhist(log=1, Plot1D=0, **kw)
        pylab.axes([v.xmin+0.40*(v.xmax-v.xmin), v.ymin+0.15*(v.ymax-v.ymin), .6*(v.xmax-v.xmin), .8*(v.ymax-v.ymin)])
        bins=10**scipy.linspace(scipy.log10(self.g_bins[0]), scipy.log10(self.g_bins[-1]), binlength) if scipy.any(scipy.diff(scipy.diff(self.g_bins))>1e-14) else scipy.linspace(self.g_bins[0], self.g_bins[-1], binlength)
        normcov = self.correlate(bins, rollavg_n=rollavg_n, normtointerval=normtointerval, scale_clim=scale_clim)
        del(pylab.gcf().axes[pylab.gcf().axes.index(parentax)])
        return normcov
    def plot_single(self, ax1, i, offset=0, PlotLog=True, **kw):
        XLIM=ax1.get_xlim()
        YLIM=ax1.get_ylim()
        ax1.plot(10*self.bjdata.delta_opening*(scipy.arange(len(self.bjdata.breake[i]))-self.breakpoints[i])+offset,
                 scipy.log10(self.bjdata.breake[i]) if PlotLog else self.bjdata.breake[i], **kw)
        ax1.set_xlim(XLIM)
        ax1.set_ylim(YLIM)
        pylab.draw()
    def makecomment(self, newpage=True):
        comments = self.bjdata.meta['Comments'] + r"\newline "+ "Speed:" + str(self.bjdata.speed) + r"\newline " + "Traces:" + str(len(self.bjdata.breake)) + r"\newline " + "Bias voltage = " + self.bjdata.meta['Bias voltage'] +r"\newline " +("normalization factor = %.2e" % self.bjdata.delta_opening) + r"\newline corrections file:" + ((self.bjdata.meta['logarithmic amplifier']).split('\\')[-1][:-2]).replace('_', '\_') +(r"\newpage" if newpage else r"\newline")
        comments = comments.replace('\r\n', '')
        return comments
        
    def doReport(self, folder, bins, basename, pars=(.08, 1e-2), xval={'breake':(-.2, 8), 'make':(-10, 0.3)}, 
                  lengthoffit=.4, rollavg_n=1, **kw):
        s=basename
        kw['cmap'] =    kw['cmap'] if kw.has_key('cmap') else pylab.cm.jet #
        self.make_subs(pars=pars, force=True)
        self.make_histograms(g_bins=bins,force=True,  breake=True)
        self.make_histograms(g_bins=bins,force=True,  breake=False)
        for bb    in ['breake', 'make']:
            f1 = pylab.figure(figsize=(11.9, 7.925));pylab.clf()
            self.plot2Dhist(breake=(bb=='breake'), xval=xval[bb], **kw)
            f1.axes[1].images[0].set_clim((-5.6, -2.2) if kw.get('log', False) else (0, 2e-3))
            mysavefig(folder, "%s_%s"%(s, bb))
            pylab.close(f1)
        f3 = pylab.figure(figsize=(11.5     ,   17.6375));pylab.clf()
        self.plot2Dhist(xval=xval['breake'], sub=True, **kw)
        pylab.figtext(0.4, 0.95, s+"_breake")
        mysavefig(folder, "%s%s_breake"%(s, '_sub'))
        pylab.close(f3)
        f11 = pylab.figure(figsize=(14, 6))
        self.plot_cycle()
        mysavefig(folder, "%s%s"%(s, '_cycle'))
        pylab.close(f11)
        f4 = pylab.figure(figsize=(9, 7.25));pylab.axes([.145, .05, .815, .815])
        self.correlate(10**scipy.linspace(scipy.log10(bins[0]), scipy.log10(bins[-1]), 200), rollavg_n=rollavg_n)
        mysavefig(folder, "%s%s_breake"%(s, '_corr'))
        pylab.close(f4)
        f4 = pylab.figure(figsize = [ 10.7,    6.        ]);pylab.axes([0, 0, 1, .9])
        self.Plot2D_corr(xval=xval['breake'], binlength=200, normtointerval=.90)
        mysavefig(folder, "%s%s_breake"%(s, 'Plot2D_corr'))
        pylab.close(f4)
        f4 = pylab.figure(figsize=(12, 7.25))
        pylab.imshow(scipy.log10(self.hist_alltraces), origin='lower', 
                     extent = [0, len(self.bjdata.breake), scipy.log10(self.g_bins[0]), scipy.log10(self.g_bins[-1])], aspect='auto')
        pylab.xlabel('Trace nr.')
        pylab.ylabel('$G/G_0$')
        pylab.grid(color='black')
        pylab.clim(0, scipy.log10(self.hist_alltraces).max()/2.)
        mysavefig(folder, "%s%s_breake"%(s, '_timeevolution'))
        pylab.close(f4)
        f5 = pylab.figure(figsize=(9, 7.25))
        self.platolength_im(plotodef_conddecade=.5, cmap=kw['cmap'] )
        mysavefig(folder, "%s%s_breake"%(s, '_plotolength'))
        pylab.close(f5)
        f6 = pylab.figure(figsize = [ 10.05,    7.25] )
        pylab.axes([0        , 0.5        , .55 , 0.5]);plotim("results/%s/%s%s.png" % (folder, s, "_breake"))
        pylab.axes([0        , 0    , .55 , 0.5]);plotim("results/%s/%s%s.png" % (folder, s, "_make"))
        pylab.axes([0.55    , 0.5    , .45 , 0.5]);plotim("results/%s/%s%s.png" % (folder, s, "_corr_breake"))
        pylab.axes([0.55    , 0    , .45 , 0.5]);plotim("results/%s/%s%s.png" % (folder, s, "_plotolength_breake"))
        mysavefig(folder, "%s%s"%('all_', s, ))
        pylab.close(f6)
        self.plot_slopes(lnm=lengthoffit, g_bins_lim = [scipy.log10(bins[0]), scipy.log10(bins[-1])], binslen=200)
        mysavefig(folder, "%s%s"%(s, 'PlotSlopes'))
        pylab.close()
        plotims([["%s/%s%s.png" % (folder, s, '_cycle'), "%s/%s%s.png" % (folder, s, '_make')]], dofigtext=False)
        mysavefig(folder, "%s%s"%('all3_', s))
        pylab.close()
        f=pylab.figure(figsize=(7.5, 6.));pylab.axes([0.125, .125, .825, .775])
        self.CalcLengthPlateaus(Glim = [.5, 1.2], bins=len(self.bjdata.breake)/20.)
        mysavefig(folder, "%s%s"%(s, '_length_of_1G0_platau'))
        pylab.close(f)
        f=pylab.figure(figsize=(7.5, 6.));pylab.axes([0.125, .125, .825, .775])
        gmax=self.g_bins[self.g_bins[:-1]<.5][self.his_breake[:, self.g_bins[:-1]<.5].sum(0).argmax()]
        self.CalcLengthPlateaus(Glim=gmax*scipy.array([.1, 10]), bins=len(self.bjdata.breake)/20.)
        mysavefig(folder, "%s%s"%(s, '_length_of_max_platau'))
        pylab.close(f)
        f4 = pylab.figure(figsize=(9, 7.25));pylab.axes([.125, .05, .835, .835])
        self.correlate_metallicwire(normtointerval=0.80)
        mysavefig(folder, "%s%s"%(s, 'correlate_metallicwire'))
        pylab.close(f4)
        f4 = pylab.figure()
        pylab.hist(scipy.concatenate(map(lambda data,firstpoint: data[firstpoint:firstpoint+10],self.bjdata.breake,self.breakpoints)), 
                   bins=10**scipy.linspace(-5, 0, 100))
        pylab.axis(xmin=.1)
        pylab.xscale('log')
        pylab.close(f4)
        f4 = pylab.figure(figsize = [ 13.7625,    7.3625])
        pylab.subplot(121);self.SnapVsPlateaulength(c_lim=(0.5, 3.2));pylab.title('0.5<G<3.2')
        pylab.subplot(122);self.SnapVsPlateaulength(c_lim=(0.5, 1.2));pylab.title('0.5<G<1.2')
        pylab.close(f4)
        pylab.close('all')
        #self.savez(self.bjdata.filenames.replace('*.dat', 'anal.dat.gz'))
            
class TUDelftmcbj(Bjdata):
    def __init__(self, **experimental):
        Bjdata.__init__(self)
        self.experimental = experimental
        for key,defaultval in [('speed', 2.),
                                ('gear',246),
                                ('v_offset', -3.e-3),
                                ('val_to_gapsize',None)]:
            self.experimental.setdefault(key,defaultval)
    def correctdata(self):
        if self.meta.has_key('Bias voltage'):
            self.experimental['Vbias'] = float(self.meta['Bias voltage'])
            correct = self.experimental['Vbias']/(self.experimental['Vbias'] + float(self.experimental['v_offset']))
            print "Data is correted to %s procent" % (correct*100)
            self.breake = map(lambda x: x*correct,self.breake)
            self.make  = map(lambda x: x*correct,self.make)
    def beforeimport(self,kw):
        """
        
        """
        if self.meta.has_key('Speed (V/s)'):
            self.experimental['speed'] = float(self.meta['Speed (V/s)'])
        if self.experimental['val_to_gapsize'] is None:
            if self.meta['logarithmic amplifier'] == 'D:\\Labview Data\\Lithographic MCBJ\\MED24_070412.TXT, ': # if new setup
                self.experimental['val_to_gapsize'] = (.28e-6*5e-5*1e9 if self.meta.has_key('Speed (V/s)') else .05)     #nm/V /motor: 0.05 nm/um
            else:
                self.experimental['val_to_gapsize'] = 0.05
        um_turn = 150.    if self.meta['logarithmic amplifier'] == 'D:\\Labview Data\\Lithographic MCBJ\\MED24_070412.TXT, ' else 250
        conv_motorcounts_to_um = (3000./um_turn*self.experimental['gear'])
        kw.update([('val_to_gapsize',self.experimental['val_to_gapsize']*self.experimental['speed']),
                   ('conv_motorcounts_to_nm',conv_motorcounts_to_um*self.experimental['val_to_gapsize'])])
        return kw
        


        
       